#include <iostream>
#include <cusp/io/matrix_market.h>
#include <cusp/monitor.h>
#include <cusp/krylov/cg.h>
#include <cusp/graph/multicoloring.h>

// where to perform the computation
typedef cusp::host_memory MemorySpace;

// which floating point type to use
typedef float ValueType;

// test function.
template< typename Matrix >
void test( )
{
    // create an empty sparse matrix structure.
    Matrix A;
    // load a matrix stored in MatrixMarket format
    cusp::io::read_matrix_market_file(A, "A.mtx");

    // Solve without preconditioning
    {
#if 0
        std::cout << "\nSolving with no preconditioner" << std::endl;
    
        // allocate storage for solution (x) and right hand side (b)
        cusp::array1d<ValueType, MemorySpace> x(A.num_rows, 0);
        cusp::array1d<ValueType, MemorySpace> b(A.num_rows, 1);

        // set stopping criteria (iteration_limit = 100, relative_tolerance = 1e-6)
        cusp::verbose_monitor<ValueType> monitor(b, 100, 1e-6);
        
        // solve
        cusp::krylov::cg(A, x, b, monitor);
#endif
    }

    // solve with diagonal preconditioner
    {
        std::cout << "\nSolving with MCSSOR preconditioner (M = )" << std::endl;

        cusp::array1d<int, MemorySpace> coloring( A.num_rows, 0 );
        cusp::graph::greedy_multicoloring( A, coloring );
        
        // allocate storage for solution (x) and right hand side (b)
        cusp::array1d<ValueType, MemorySpace> x(A.num_rows, 0);
        cusp::array1d<ValueType, MemorySpace> b(A.num_rows, 1);

        // set stopping criteria (iteration_limit = 100, relative_tolerance = 1e-6)
        cusp::verbose_monitor<ValueType> monitor(b, 100, 1e-6);
        // setup preconditioner
        // cusp::precond::diagonal<ValueType, MemorySpace> M(A);
        // solve
        // cusp::krylov::cg(A, x, b, monitor, M);
    }
}

int main(void)
{
    // create an empty sparse matrix structure.
    typedef cusp::coo_matrix<int, ValueType, MemorySpace> Matrix;
    test<Matrix>( );
    return 0;
}

